﻿ ESTIMAREA CU PLAUZIBILITATE MAXIMĂ, ML În unele situaţii fie nu putem să găsim un estimator MVU, fie el nu există În astfel de situaţii putem să recurgem la estimarea cu plauzibilitate (verosimilitate) maximă Estimatorii ce rezultă dintr-o astfel de abordare sunt cei mai răspândiţi în aplicaţiile practice, în rezolvarea unor probleme complicate de estimare În multe cazuri, performanţele acestor estimatori sunt aproape optimale, dacă înregistrările de date sunt foarte lungi În consecinţă, estimatorii ce maximizează plauzibilitatea, sunt aproximări bune ale estimatorilor MVU Dacă estimatorul este suboptimal dar ˆ limEθ θ= {} N→∞ 1− ˆ limDisp CRLB Iθ θ== () {} Nθ →∞ se spune că estimatorul este asimptotic nedeplasat, şi asimptotic eficient Pentru înregistrări lungi, adică N destul de mare, un astfel de estimator poate fi considerat ca fiind optimal 1 Determinarea estimatorului de plauzibilitate maximă Estimatorul de plauzibilitate maximă, ML, este acea valoare a parametrului necunoscut, θ, care maximizează p(x; θ), cu x fixat, ceeace înseamnă că maximizează funcția de plauzibilitate, L(θ) ˆ ==x=x arg max max ;Lpθ () () MLθθ 0 θθ Probabilitatea ca x observat, să cadă într- un volum Δx ce tinde spre 0, centrat pe x0 este p(x=x0; θ)Δx şi ea depinde, în principiu, de valoarea adevărată a parametrului θ În figură se arată p(x=x0; θ) ca funcţie de θ Dacă valoarea adevărată a parametrului θ ar fi θ 3atunci probabilitatea de a se fi observat x=x0ar fi 3 este una puțin credibilă, foarte mică, practic nulă Prin urmare, ipoteza θ=θ puțin plauzibilă Ipotezele θ=θ 2 sau θ=θ4, în condițiile în care s-a observat că x=x0,vor fi respinse pe aceleaș i motive Urmând acelașiraționament, ipoteza 1 este cea mai plauzibilă, deoarece p(x=x0; θ=θ1) are valoarea maximă θ=θ Pentru exemplificare vom considera un semnal cu modelul xnAwn A wn An N xn AA=+ > = −∼∼NN; 0, 0, 0,1, , 1 , ()() [][][][] i dispersia zgomotului alb, fiind atât valoarea componentei continue cât ș2A gaussian 1 Densitatea de probabilitate a datelor are expresia 1N− 211⎧⎫ x −> () [] () ⎨⎬ ∑ /2; exp ; 0NpA xnA A=− 022nAAπ = ⎩⎭ () Dacă datele sunt fixate, adică cele N componente ale vectorului x se înlocuiesc cu cele N valori efectiv măsurate, se obține funcția de plauzibilitate, o funcție ce depinde doar de parametrul necunoscut, A În formă logaritmică funcția este 1N− 21NN 0lA p A A xn A Aπ==−−− − >xxfixatln ; ln 2 ln ;, () ( )() [] () ∑ 0222nA= Maximul ei apare pentru A>0, pentru care derivata este nulă 11NN−− 211dl AN () + − + − = 0xn A xn A=− [][] ()() 2∑∑ 0022nndA A AA== Dezvoltând pătratul și efectuînd calculele, se obține ecuația 11 1NN N−− − 211 1NN − ++ −= 0xn xn xn N−+ [][][] ∑ 2∑∑ 00 0222nn nAAAA== = 3 sau 1N− 221 = > 0; 0AA xn A+− [] ∑ 0nN= Soluțiile ecuației sunt 1N− 211 1 ± + [] ∑ 1,2Axn=− 024nN= dar numai soluția pozitivă are sens în exemplu 1N− 211 1 ˆ Axn=−+ + [] ML∑ 024nN= 4 2 Vom considera modelul de semnal în care nu cunoaştem componenta continuă A, zgomotul w[n] fiind alb, gaussian cu dispersia cunoscută 2 ; 0,1, , 1; 0,xn A wn n N wnσ=+ = −∼N [][][] () Densitatea de repartiție a datelor x este 1N− 211⎧ ⎫ x − () [] () ⎨⎬ ∑ 2/2;expNpA xnA=− 22nσ 0= ⎩⎭ 2πσ () Cu x fixat (adică înlocuit cu vectorul obținut prin măsurare) plauzibilitatea logaritmică, dependentă doar de parametrul necunoscut, A, este 1N− 21N 2 ln ; ln 2 ; fixatlA p A xn Aπσ ==− − −xx () ( ) [] () () 2∑ 022nσ = Ecuația ce permite daterminarea estimatului ML pentru A se obține anulând derivata, în raport cu A, a plauzibilității logaritmice 11NN−− () 11lA∂ =−= 0xn A xn NA=− [][] () 22∑∑ ∂ 00nnAσσ == Se obține estimatul ML, care în acest caz, este chiar media eșantion 1N− 1 ˆ = [] ∑ MLAxn x= 0nN= 5Dacă într-o problemă există un estimator eficient, atunci procedura ML îl va produce ! Proprietăți asimptotice ale estimatorilor ML Estimatorii ML sunt, asimptotic, nedeplasați șieficienți, adică 1ˆˆ− lim ; limEDispCRLBIθ θθ θ=== () {}{} →∞ →∞ NNθ Simbolic, distribuția asimptotică a estimatorului ML se notează cu a 1− ˆ θθN ,Iθ () () ∼ Pentru primul exemplu, în care atât media cât și dispersia datelor este A, primele două derivate ale plauzibilității logaritmice sunt 11NN−− x 2ln ;11pAN∂ () + − + − [][] ()() ∑∑ 2xn A xn A=− ∂ 0022nnAAAA== 211 1NN N−− ∂x 2ln ;11 1pAN − () − −− − [][][] ()() ∑∑ ∑ 222 2 3xnxnAxnA=− ∂ 00 02nn nAAA A A== = Mediind statistic a doua derivată obținem 21N− ln ;1pAN⎧⎫x ∂ () ⎪⎪ [] {} ⎨⎬ ∑ 222EExn=− ∂ 02nAAA= ⎪⎪ ⎩⎭ 11NN−− 211 −− [][] {}() () ∑∑ {} 23Exn A E xn A−− 00nnAA== NNN NN⎛⎞ =− + =−6 ⎜⎟ 222− 22AAAAA⎝⎠ 3 În consecință inversa informației Fisher este 2 1111A− IA=− = = () 2 1NNN⎧⎫ x ln ;pA∂ () A⎪⎪ ++ 2E⎨⎬ 222AA A∂ ⎪⎪ ⎩⎭ Am stabilit deci repartiția asimptotică a estimatorului ML pentru A 2a 1A⎛⎞ ˆ ,AAN ⎜⎟ ∼ ML 12NA+ ⎝⎠ Pentru A=1 avem, 112− NI A= = () 1112 3A= + produs constant, independent de N Pentru a ne forma o idee despre ce înseamnă “N suficient de mare” am realizat un experiment pentru lungimi tot mai mari ale datelor, N=5, 10, 15, ,50 Am repetat de M=10000 de ori experimentul, calculând 10000 de valori ale mediei și ale dispersiei Pentru a reduce fluctuațiile statistice am calculat medii ale mediilor și ale dispersiilor obținând estimări ale medie statistice și ale dispersiei (statistice) 2 MM⎛⎞ 11^^ ^ ˆˆˆ ˆˆ −= EA A DispA A EA M⎜⎟; ; 10 000 {}{} ∑∑ ii== 7{} ⎜⎟ 11iiMM== ⎝⎠ Se observă că estimarea mediei ne Produsul N Dispersia este practic la sugerează că pentru N>25 ne-am valoarea 2/3 pentru N>25 apropiat suficient de valoarea adevărată A=1 pentru a declara că estimatorul ML este nedeplasat În mod uzual, N suficient de mare pentru a admite că este adevărată repartiția , înseamnă câteva zeci de eșantioane 8asimptotică 4 Putem da, fără demonstrație, teorema Dacă densitatea de probabilitate a datelor x, p(x; θ), satisface condiția de regularitate ;pθ∂ x ⎧⎫ () 0E= ⎨⎬ θ ∂ ⎩⎭ atunci estimatorul pentru parametrul necunoscutθ, ce maximizează plauzibilitatea, adică estimatorul ML, este distribuit asimptotic conform cu a 1− ˆ ,Iθ θθN () () ∼ I(θ) fiind informația Fisher, evaluată la valoarea adevărată a parametrului θ *** Se poate deci spune că, în mod asimptotic, estimatorul ML este nedeplasat şi eficient, adică el este asimptotic optimal Viteza de convergenţă, adică valoarea N după care se pot aplica formulele asimptotice, este o problemă deschisă În majoritatea cazurilor N nu este foarte mare, fiind, în mod uzual, de ordinul 9zecilor Estimator ML pentru faza iniţială a unei sinusoide Reluăm problema determinării fazei inițiale pentru o sinusoidă, afectată de un zgomot alb, gaussian Modelul datelor este 2 =+Φ+ =−∼N [][][] () () 0cos 2 ; 0, ; 0,1, , 1xn A fn wn wn n Nπσ Dacă se fixează x în expresia densității de probabilitate a datelor, se obține plauzibilitatea, funcție numai de parametrul necunoscut, Φ 1N− 211⎧⎫ ⎡⎤ Φ= − − +Φ x () [] () ⎨⎬ ∑ 02/2exp cos 2 ; fixatNLxnAfnπ ⎣⎦ 22σ 0n= ⎩⎭ 2πσ () Pentru stabilirea estimatorului ML e necesar să maximizăm L(Φ), ceeace e echivalent cu a minimiza, după Φ, funcția 1N− 2 ⎡⎤ Φ= − +Φ cos 2JxnAfnπ () [] () ∑ 0⎣ ⎦ 0n= Anulăm derivata acestei funcții 1N− dJΦ () ˆˆ⎡⎤ 2cos2 sin20xn A fn A fnππ = − +Φ +Φ = [] ()() 00∑ 10d⎣⎦ Φ 0n= 5 Rezolvarea ecuației, aplicând aproximările cunoscute, conduce la 11NN−− A ˆˆ +Φ= +Φ≅ [] ()() ∑∑ 00sin 2 sin 4 2 0xn fn fnππ 002nn== cu condiția ca frecvența digitală să nu fie foarte apropiată de 0 sau 0 5 Dacă dezvoltăm membrul stâng al relației de mai sus obținem 1N− ˆ +Φ = [] () ∑ 0sin 2xn fnπ 0n= 11NN−− ˆˆ Φ+Φ≅ [] [] ∑∑ 00sin 2 cos cos 2 sin 0xn fn xn fnππ 00nn== din care se determină estimatorul ML pentru faza inițială ca fiind 1N− [] ∑ 0sin 2xn fnπ x () 20ˆnT= − =− arctg arctgΦ≅ 1MLN− () 1Tx [] ∑ 0cos 2xn fnπ 0n= 11 relație în care am utilizat notațiile cunoscute deja 11NN−− ==xx ()() [][] ∑∑ 1020cos 2 ; sin 2TxnfnT xnfnππ 00nn== În Cap 2 am determinat informația Fisher pentru faza inițială a sinusoidei 2 NA IΦ≅ () 2 2σ Aplicând teorema enunțată putem afirma că repartiţia asimptotică a estimatorului fazei este normală: 2a ⎛⎞ 2σ ˆ N ,ΦΦ ∼ ML⎜⎟ 2 NA ⎝⎠ Dispersia estimatorului ML al fazei tinde asimptotic spre 22 21Aσ ˆ Φ→ = == Disp SNRη ; {} 22ML σ 2NNAη 12 6 Dacă introducem expresia semnalului de date în relația estimatorului ML al fazei inițiale obținem 1N− ⎡⎤ +Φ + [] () ∑ 00cos 2 sin 2Afnwn fnππ ⎣⎦ 0ˆn= Φ≅− 1MLNarctg− ⎡⎤ +Φ + [] () ∑ 00sin 2 cos 2Afnwn fnππ ⎣⎦ 0n= 11 1NN N−− − −Φ+ +Φ+ [] () ∑∑ ∑ 00sin sin 4 2 2 sin 2AAfn wnfnππ 00 0nn n== = ≅− 11 1NN Narctg−− − Φ+ + Φ + [] () ∑∑ ∑ 00cos cos 4 2 2 cos 2AAfn wnfnππ 00 0nn n== = Dacă frecvența digitală nu este foarte aproape de 0 sau 0 5, termenii al doilea de la numărător și de la numitor pot fi neglijați, așa că rezultă 1N− −Φ+ [] ∑ 0sin 2 sin 2NAwnfnπ 0ˆn= Φ≅− 1MLNarctg− Φ+ [] ∑ 0cos 2 cos 2NAwnfnπ 0n= 1N− 2 Φ− [] ∑ 0sin sin 2wn fnπ 0nNA= = Narctg−1 2 Φ+ ∑ 13[] 0cos cos 2wn fnπ 0nNA= Estimarea cu plauzibilitate maximă în cazul transformării parametrilor În unele cazuri nu ne interesează valoarea parametrului necunoscut, θ, de care depinde în mod nemijlocit semnalul de date, ci un parametru derivat, α= g(θ) De exemplu, pentru datele, în care parametrul necunoscut este A 2 ; 0,1, , 1; 0,xn A wn n N wnσ=+ = −∼N [] [][] () dorim să estimăm A exn wnαα==+; ln [][] Cu datele x fixate, funcția de plauzibilitate dependentă de parametrul A este 1N− 211⎧⎫ x −∈ () [] () ⎨⎬ ∑ NpA xnA A=−2/2;exp ; 22σ 0n= ⎩⎭ 2πσ () Această funcție de plauzibilitate se poate se poate transforma ușor într-o funcție de plauzibilitate dependentă de parametrul α, deoarece, în acest caz funcția g( ) este biunivocă 1N− ⎫ 211⎧ =−− x ;exp lnpxnαα () [] () ⎨⎬ ∑ α 142/22N 02nσ = ⎩⎭ 2πσ () 7 Se maximizează, în raport cu α, această funcție dependentă de α 1N− ln ;11dpαx () α ˆlnxnα =− [] () 2∑ σ = 0ˆndαα 1N− 1⎛⎞ ln 0xn Nαˆ =−= [] ∑ 2⎜⎟ = 0ˆnασ ⎝⎠ Din anularea parantezei rezultă logaritmul estimatorului ML pentru parametrul α 1N− 1 ˆlnxnxα == [] ML∑ 0nN= din care se obține estimatorul ML pentru parametrul α, obținut prin transformare x ˆeα = ML Se știe că media eșantion este estimatorul ML pentru componenta continuă, așa că putem scrie că ˆ ˆMLA ˆegAα () MLML== 15Relaț ia de transformare a parametrilor este șirelația de transformare a estimatorilor ML, cel puțin pentru transformările g( ) biunivoce Vom da următoarea teoremă fără demonstrație Estimatorul ML pentru gαθ= () unde θ parametrizează densitatea de probabilitate a datelor, p(x;θ), este dat de relația ˆ ˆgαθ () MLML= ˆ Estimatorul se obține maximizând după θfuncția de plauzibilitate p(x;θ), cu MLθ x fixat *** Vom considera un exemplu în care datele sunt un zgomot alb, gaussian, cu puterea necunoscută 2 ; 0,1, , 1; 0,xn wn n N wnσ==−∼N [][][] () pentru care dorim să estimăm puterea în unități logaritmice 2 dBPσ1610log = [] 8 Deoarece zgomotul este alb, cu x fixat, obținem funcția de plauzibilitate 1N− ⎫ 2211⎧ pxnσ;exp =− x [] ()⎨⎬ 2/2N∑ 22σ n= 0⎩⎭ 2πσ () Se determină estimatorul ML pentru putere anulând derivata formei logaritmice 2 − 1ln ;Ndpσ x () 21N =0xn=− + [] 224∑ σ = 022ndσσ Din care se determină estimatorul ML pentru puterea zgomotului 1^N− 221 = [] ∑ MLxnσ 0nN= Substituim rezultatul în relația de transformare șiobținem noul estimator ML 1N− 21⎛⎞ ˆ Pxn dBα10 log 10 log () [] [ ] ∑ ⎜⎟ ML ML== 0nN= ⎝⎠ 17 Determinarea numerică a estimării ML Un avantaj al estimării cu plauzibilitate maximă este acela că putem determina o valoare estimată pentru fiecare set de date În fapt, se caută maximul unei funcţii (de plauzibilitate), complet definită de un set de date Dacă ştim că parametrulθ aparţine unui interval precizat, aşa cum se arată în figură, putem recurge la o împărţire destul de fină a intervalului, printr-o grilă Calculând valorile plauzibilităţii în punctele alese, putem stabili un argument ce maximizează , plauzibilitatea Putem să recurgemşi la subîmpărţirea intervalului în care am găsit un maxim, în aşa fel încât localizarea punctului A (de maxim) să fie destul de bună Există riscul să găsim un punct de maxim local, cum este B, în loc de maximul absolut, A Este evident că pentru date noi, procedura de calcul trebuie reluată, deoarece valoarea estimatorului ML se 18schimbă 9 Dacă domeniul valorilor posibile ale parametruluiθnu este un interval finit, căutarea maximului prin grilă pune serioase probleme de calcul În astfel de cazuri se apelează la metode iterative de determinare a valorii argumentului ce maximizează plauzibilitatea Dintre metodele utilizate vom prezenta, pe scurt, algoritmul Newton-Raphson Modificarea datelor impune reluarea calculului, de fiecare dată Trebuie specificat faptul că alegerea valorii iniţiale influenţează mult convergenţa algoritmului Pentru exemplificare, vom considera modelul de date x[n], în care zgomotul w[n] este alb, gaussian iar parametrul necunoscut, r, este pozitiv 2n ; 0,1, , 1; 0,xn r wn n N wnσ=+ = −∼N [][][] () Dacă în expresia densităţii de probabilitate (repartiţie) fixăm vectorul de date x, obţinem funcţia de plauzibilitate, dependentă doar de parametrul necunoscut, r 1N− 211⎧ n⎫ x pr xnr=−;exp − () [] () ⎬ ∑ 2/2N⎨ 22σ n= 0⎩⎭ 2πσ () 19 Căutăm valoarea parametrului r care maximizează plauzibilitatea Maximizarea plauzibilităţii este echivalentă cu minimizarea, după r, a funcţiei 1N− 2 n Jrxnr=− () [] () ∑ 0n= Estimatorul ML este soluţia ecuaţiei 1N− () 1nndJ r− ˆˆ20xn r nr=− = [] () ∑ 0ndr= O astfel de ecuaţie nu poate fi rezolvată, în sensul stabilirii unei formule se calcul a soluţiei *** Vom prezenta metode iterative de maximizare a plauzibilităţii logaritmice, ce caută un zerou al derivatei sale (o rădăcină): x ln ;pθ∂ () 0= ∂ θ Soluţia ecuaţiei se determină în mod iterativ (prin mai mulţi paşi), pornind de la funcţia o valoare iniţială care se adoptă Se notează cu g(θ) x ln ;dpθ () = gθ 20() dθ 10 Algoritmul Newton-Raphson Considerăm că valoarea de pornire în iteraţia pentru căutarea rădăcinii este θ0 Ecuaţia tangentei la curba y=g(θ ),dusă în punctul B (vezi figura) este dgθ () θθ−= − ()() 00ygθ dθ = 0θθ Intersecţia acestei tangente cu axa Ox, ce se obţine punând y=0, defineşte o soluţie mai bună pentru ecuaţia ce trebuie rezolvată Ea este θ () () 0gdgθ =− = θ 21() 10 0; 'gθθ θθ= θ () 00'dgθ Se duce tangenta în punctul D şi se obţine o nouă soluţie, corespunzătoare punctului E gθ () 1 θθ=− 21 'gθ () 1 se determină cu relaţia de iterare Valorile succesive ale soluţiilor pentru θ gdgθ () () kθ ; 'gkθθ=− = = 0,1, 2, θ () 1kk k+ = θ 'dgθθ () kkθ Condiţia de oprire a procesului iterativ de determinare a soluţiei este θθε−< 1kk+ care afirmă că soluţiile succesive se află într-un interval cu lungimea mai mică decât eroarea maximă admisă ε Substituind forma funcţiei g(θ)în formula de iterare, algoritmul Newton- Raphson pentru determinarea iterativă a estimatorului ML pentru θ devine 1− 2 ⎡⎤ ln ; ln ;dp dpθθ xx ()() =− = 0,1, 2, kθθ ⎢⎥ 12kk+ = θθ θ ⎢⎥ kddθ 22⎣⎦ 11 Vom face câteva precizări privind aplicarea algoritmului Newton- Raphson: i) Este posibil ca şirul de valori obţinut prin iteraţie să nu conveargă Acest lucru se întâmplă mai ales atunci când derivata a doua a plauzibilităţii logaritmice este de valoare apropiată de zero; drept consecinţă, termenul de corecţie din relaţia (6 61) fluctuează mult, de la iteraţie la iteraţie ii) Chiar dacăşirul de valori obţinut prin iteraţie converge, s-ar putea să nu ne conducă la un maxim global (punctul A din figură) ci la un maxim local (punctul B din figură), fie chiar într- un punct de minim (punctul C din figură) Pentru a evita cantonarea într-un maxim local este bine să lucrămcumai multe valori iniţiale şi să reţinem acea valoare care conduce la “cel mai mare maxim” al funcţiei de plauzibilitate Alegerea valorii iniţiale este deosebit de importantă într-o astfel de procedură iterativă 23 Revenim la exemplul cu care am început Avem, pentru derivatele plauzibilităţii logaritmice 1N− x () 1ln ;1nndp r− xn r nr=− [] () 2∑ 0ndrσ = 21N− x () 2ln ;1nndpr− ⎤ =−−− 121nr n x n n r⎡ ()() [] ∑ 22⎣ ⎦ 0ndrσ = Iteraţia pentru stabilirea estimatului ML pentru parametrul necunoscut, r, devine 1N− 1nn− xn r nr− [] () kk∑ 0n= rr=− 11kkN+ − 2nn− ⎤ 121nr n x n n r⎡ −−− ()() [] ∑ kk⎣ ⎦ 0n= 24 12 O versiune a metodei Newton-Raphson face uz de legea numerelor mari, care afirmă că pentru eşantioane IID media aritmetică este o aproximare bună a mediei statistice, dacă numărul de eşantioane mediate este mare 25 Conform legii numerelor mari, dacă N este suficient de mare, putem scrie, pentru cazul a N eşantione de date x[n], statistic independenteşi identic distribuite (IID) că 22 1N− ln ;ln ;dpxndpθ x θ [] () () ∑ 22= 0nddθθ = 2 1N− ln ;1dpxnθ [] () ∑ N=2 0nNdθ = ↓ ↓ 2 ⎧⎫ x ln ;dpθ () ⎪⎪ ⎨⎬ NE≅2 dθ ⎪⎪ ⎩⎭ =− =− Ni Iθθ () () Pentru cazul în care N este suficient de mare şi eşantioanele de date sunt IID, formula de iterare a algoritmului Newton-Raphson devine x () 1ln ;dpθ− ; 0,1, 2, Ikθθ θ =+ = () 1kk k+ θθ = kdθ prin această relaţie este cunoscută şi sub denumirea de 26Metoda descrisă “Method of Scoring” 13 Revenim (a doua oară) la exemplul de început Pentru cazul în speţă, informaţia Fisher este 2 ln ;dpr⎧⎫x () ⎪⎪ IEθ =− () 2⎨⎬ dr ⎪⎪ ⎩⎭ 1N− 21nn− ⎤ =− − − − 1 2 1nr n E x n n r⎡ ()() [] {} ∑ 2⎣ ⎦ 0nσ = 1N− 21nnn− 1 2 1nr nr nr⎡⎤ =− − − − ()( ) ∑ 2⎣⎦ 0nσ = 1N− 22 21n− nr= 2∑ 0nσ = Dacă o substituim în formula de iteraţie “method of scoring” obţinem 1N− 1nn− xnrnr− [] () kk∑ 0n= rr=+ 11kkN+ − 22 2n− 27nr∑ k 0n= Pentru cazul valorii reale r=0 5 şi pentru dispersia de 0 01, s-a efectuat o simulare pentru estimarea MLE cu relaţia ce implementeză “method of scoring”, pentru o realizare a datelor cu lungimea N=50 Condiţia de oprire a simulării a =0 0001 În figură se pot urmări valorile r succesive, dacă valoarea iniţială fostε a fost 0 9 După numai 6 iteraţii se constată că valorile succesive au primele patru zecimale egale, aşa că valoarea estimatului ML este 0 4969 Dacă valoarea iniţială e luată 0 3 se constată că după numai patru iteraţii se obţine condiţia de oprire Valoarea estimării este aceeaşi deoarece datele sunt aceleaşi 4− 4− −< −< 110kkrr+ 110kkrr+ 28 14 Dacă însă valoarea iniţială este 1 2, condiţia de oprire se atinge doar după 17 iteraţii 4− −< 110kkrr+ Se poate observa că valoarea iniţială influenţează mult numărul de iteraţii necesare pentru calculul estimării ML 29 O extindere pentru cazul parametrului vector Reluăm problema determinarii a doi parametri necunoscuţi, dacă modelul de semnal corespunde unei componente continue, afectată de un zgomot alb, gaussian 2 ; 0,1, , 1; 0,xn A wn n N wnσ=+ = −∼N [][][] () Vectorul parametrilor necunoscuţi este T 2 ⎡⎤ θ = Aσ ⎣⎦ Cu datele x fixate, se determină plauzibilitatea Pentru a stabili estimatorul ML, se anulează derivata plauzibilităţii logaritmice, în raport cu A 1N− ln ;1p∂xθ () 0xn A= −= [] () 2∑ ∂ = 0nAσ şi se obţine 1N− 1 ˆ Axn x== [] ML∑ 0nN= Se anulează apoi derivata în raport cu al doilea parametru 1N− xθ 2ln ;1pN∂ () 30xnA=− + − [] () 224∑ ∂ σ = 022nσσ 15 Conform celor prezentate anterior, valoarea pentru A se înlocuieşte cu estimarea sa ML ˆ AAx=← ML şi se obţine ecuaţia 1N− 21N 0xn x−+ −= [] () 24∑ 022nσσ = Găsim al doilea estimator ML de forma 1^N− 21 2 xnxσ=− [] () ML∑ 0nN= Cei doi estimatori scalari se pot grupa într-un estimator ML vector x⎡ ⎤ ⎢⎥ 1ˆN− θ 21= ML ⎢⎥ xn x− [] () ∑ 31⎢ ⎥ 0nN= ⎣⎦ Dăm, fără demonstraţie, teorema ), ce depinde de parametrul Dacă densitatea de probabilitate a datelor, x, p(x; θ necunoscut, satisface condiţia de regularitate vector θ, ;p∂ xθ ⎧⎫ () E= 0 ⎨⎬ ∂θ ⎩⎭ atunci estimatorul ML pentru parametrul vector necunoscut, θ, ce maximizează plauzibilitatea, este repartizat asimptotic, conform cu a 1− ˆ ,θθ N () () MLI θ ∼ unde I(θ ) este matricea de informaţie Fisher, cu elementele evaluate la valorile adevărate ale componentelor vectorului θ *** 32 16 Reluăm problema determinarii a doi parametri necunoscuţi, dacă modelul de semnal corespunde unei componente continue, afectată de un zgomot alb, gaussian T 22 ⎡⎤ θ∼ xn A wn n N wn AσσN; 0,1, , 1; 0, ; =+ = − = [] [] [] () ⎣⎦ pentru care tocmai am stabilit estimatorii ML Cunoaştem repartiţia mediei 2 ⎛⎞ σ ,Ax∼N ⎜⎟ N ⎝⎠ Mai ştim că statistica T are o repartiţie hi-pătrat, şi este statistic independentă de media eşantion 2 1N− ⎛⎞ [] 2xn x− ∼ = ∑ ⎜⎟ 1NTχ − 0nσ = ⎝⎠ Literatura afirmă că repartiţia hi-pătrat converge spre una normală 2 ,2NNχN () N→ N→∞ Avem deci pentru statistica T o repartiţie asimptotică normală a 332 1, 2 1TNNχ−−∼∼N () () 1N− Forma estimatorului ML pentru dispersia datelor a fost stabilit de noișieste 2 2211^NN−− − 21xn xσσ [] 2⎛⎞ =−= = [] () ⎜⎟ ∑∑ MLxnx Tσ 00nnNN Nσ == ⎝⎠ Media estimatorului se determină cu relația 22^⎧⎫ ⎪⎪ 2σσ 1EETNσ==− {}() ML⎨⎬ NN ⎪⎪ ⎩⎭ iar dispersia are forma 44^⎧⎫ ⎪⎪ 2σσ 21Disp Disp T Nσ==− {}() 22ML⎨⎬ NN ⎪⎪ ⎩⎭ Se poate da acum repartiția asimptotică a estimatorului ML pentru dispersie ^ a ⎛⎞ 22411NN−− σσ∼N ,2σ 2ML⎜⎟ NN⎝⎠ Cele două componente ale vectorului estimator sunt repartizate normal șisunt statistic independente Componentele vectorului sunt deci și mutual gaussiene Caracterizarea sa statistică se realizează, în consecință, doar prin vectorul 34 mediilor și prin matricea de covarianță 17 Vectorul mediilor demonstrează o comportare asimptotică de neabatere Ex⎡⎤ {} A⎡⎤ A⎢⎥ ⎡⎤ ˆ^⎧⎫ ⎢⎥ EN⎢⎥ θθ = {} ⎢⎥ ⎪⎪ → ML==222 ⎢⎥ ⎢⎥ σ →∞ NEσ ⎨⎬ ⎣⎦ MLσ 1N− ⎣⎦ ⎢⎥ ⎪⎪ ⎩⎭ ⎣⎦ Matricea de covarianță are pe diagonala principală dispersiile Celelalte două elemente sunt nule, ca urmare a lipsei corelației între cei doi estimatori ML scalari 22 ⎡⎤ ⎡⎤ σ σ 00 ⎢⎥ ⎢⎥ 1NN− ⎢⎥ CθI θ == ⎢⎥ ()() 4→ − ⎢⎥ →∞ ⎢⎥ 412NNσ 020σ 2⎢⎥ NN⎢⎥ ⎣⎦ ⎣⎦ Se vede imediat că matricea de covarianță tinde asimptotic spre inversa matricei Fisher, stabilită de noi anterior Am verificat, prin calcul direct, faptul că teorema teorema privind repartiția asimptotică a estimatorului vector ML se aplică în exemplul de faţă: a 1− ˆ ∼N ,θθ () 35() MLI θ Teorema de invarianță a estimatorului ML vector Estimatorul de plauzibilitate maximă (ML) pentru vectorul α 1r× obținut din vectorul θ 1p× g(θ) este prin funcția inversabilă α= ˆ ˆg= () MLMLαθ Dacă funcțiaα =g(θ) nu este inversabilă, atunci estimatorul ML pentru α asigură maximizarea funcției de plauzibilitate transformată maxpp=x;αx;θ ()() T θα=gθ: () {} *** 36 18 Vom exemplifica cu modelul de semnal componentă continuă afectată de un zgomot alb, gaussian, pentru care nu se cunosc A și puterea zgomotului 2 ; 0,1, , 1; 0,xn A wn n N wnσ=+ = −∼N [][][] () dar ne interesează doar raportul semnal/zgomot, SNR 2 A α= 2 σ Având în vedere cele expuse avem pentru estimatorul SNR 2 2ˆ () () MLAx ˆ= =α 1MLN− 21^ 2xnx− [] () ∑ σ 0MLnN= Vom aborda acum cazul modelului gaussian generalizat, în care zgomotul este ) gaussian, avand media nulă şi matricea de covarianţă C(θ ,w0Cθ∼N () () Repartiția datelor x este normală 37 ,xμθCθ∼N ()() () Se poate arăta că pentru repartiţia normală de mai înainte, este valabilă regula de derivare a plauzibilității T ⎧⎫ x;θCθμθ ∂∂∂ ()()() ⎪⎪ 111p−− C θC θ x μθ tr=− + − ⎡⎤ ()() () ⎨⎬ ⎣⎦ 2θθ ∂∂∂ ⎪⎪ ⎩⎭ kkkθ 1− C θ () 1T∂ ⎡⎤ ; 1, 2, ,kp⎡⎤ −− − = x μθx μθ ()() ⎣⎦ ⎣⎦ 2θ k∂ Modelul liniar generalizat este descris de relaţia: x=Hθ +w HNxpfiind matricea de observare După ce datele x au fost fixate, expresia funcției de plauzibilitate este ⎧⎫ 111T− x;θx-Hθ Cx-Hθ expp=− ()()() ⎬ /21/ 2N⎨ 2⎩⎭ C 2π () Maximizarea plauzibilităţii implică minimizarea funcţiei: 1T− J=θx-Hθ Cx-Hθ ()()() 38 19 Putem să recurgem direct la maximizarea plauzibilităţii, ţinând seama de regula de derivare anterior expusă Avem, pentru medie =μθHθ () șideci T ∂∂x;θHθ ()() 1lnp− ; 1, 2, ,kp==Cx-Hθ () ∂∂ kkθθ Dacă scriem în mod explicit cele p relații, acestea pot fi puse sub formă matriceală T ⎡⎤ Hθ () ⎡⎤ x;θ ∂ lnp∂ () ⎢⎥ ⎢⎥ ∂ ∂ θ ⎢⎥ 11θ ⎢⎥ T⎢⎥ x;θ ∂ lnp⎢⎥ () Hθ ∂ () ⎢⎥ ⎢⎥ 1− ⎢⎥ ∂ Cx-Hθ = ⎢⎥ () ∂ θ 22θ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ⎢⎥ ∂ x;θ lnTp⎢⎥ () ⎢⎥ ∂ Hθ () ⎢⎥ ⎢⎥ ∂ ⎢⎥ 39pθ ∂ ⎣⎦ ⎢⎥ pθ ⎣⎦ sau, cu notațiile matriceale uzuale T ∂∂x;θHθ ()() 1lnp− =Cx-Hθ () ∂∂θθ Cunoaștem regula de derivare T Hθ () T∂ = H ∂ θ care, după aplicare, ne dă ∂x;θ () 1lnTp− =HC x-Hθ () θ ∂ Egalăm derivata cu 0 șiobținem estimatorul ML pentru vectorul θ 11TT−− ˆ HC HθHC x ML= cu expresia 1− 11TT− ˆ− θHC H HC x () ML= Rezultă că, pentru modelul liniar generalizat, estimatorul de plauzibilitate maximă esteşi estimator MVU şi este şi eficient ˆˆ =θθ 40MLMVU 20 Estimatorul are matricea de covarianță de forma 1− 1T− =CHCH () ˆ θ 41 Vom considera că într-o problemă de estimare, în care datele x conduc la ), există un estimator eficient funcția de plauzibilitate p(x; θ ˆ efθ Fiind eficient, el satisface teorema Cramer-Rao, adică x;θ lnp∂ () ˆ −Iθθ θ () () ef= ∂θ Estimatorul ML satisface condiția de maxim a plauzibilității logaritmice x;θ lnp∂ () 0 = ˆ θ θθ = ML∂ Substituind valoarea estimatorului ML în expresia teoremei Cramer-Rao obținem x;θ lnp∂ () ˆˆˆ =−=Iθθθ0 () () ˆML MLef θ =θθ ML∂ Rezultă concluzia ˆˆ θθ MLef= Concluzia se poate enunța ca teoremă Dacă într-o problemă de estimare există un estimator eficient, atunci procedura 42 de determinare a estimatorului ML îl va produce, cei doi estimatori fiind egali 21 Estimarea asimptotică cu plauzibilitate maximă În literatura de specialitate se arată că, pentru N suficient de mare, adică asimptotic, plauzibilitatea logaritmică se poate aproxima cu relația 1/ 2 ⎤ IfNN⎡ () ln ln 2 lnpPfdfπx;θ − + ()() ⎢⎥ xx≅− ∫ () ⎢⎥ 1/ 222xxPf− ⎣⎦ relație în care Pxx(f) este densitatea spectrală de putere a datelor iar I(f) este periodograma datelor x 2 1N− 211 − 2jfnπ Ifxne Xf== ()() [] ∑ 0nNN= are loc prin Pxx(f) Derivata Dependența plauzibilității de parametrul vector θ plauzibilității logaritmice după componenta θ ise calculează cu relația 1/2 ln1pIfPfN⎡⎤x;θ ()()() xx∂∂ − = ; 1, 2, ,df i p≅− ⎢⎥ ∫ ∂∂ ()() 1/22xx xxiiPf Pfθθ − ⎣⎦ Pentru stabilirea celor p componente ale estimatorui vector ML, este necesară rezolvarea sistemului de ecuații 1/2 ∂ ()() 1xxIf P f⎡⎤ 0; 1, 2, ,dfip−= = ∫ 43⎢⎥ ∂ ()() 1/2xx xxiPf Pfθ − ⎣⎦ Câteva exemple din domeniul prelucrării semnalelor Estimarea distanţei în RADAR şi în SONAR Am determinat, pentru RADAR, densitatea de repartiție a datelor Dacă în expresia acesteia se consideră datele fixate, se obține funcția de plauzibilitate 012n− 1xn⎧⎫ [] ⎪⎪ x () ⎨⎬ ∏ 02;exp pn=− πσ = 022nσ ⎪⎪ ⎩⎭ 2 ⎧⎫ 01nMxn sn n+− [] [ ] () ⎪⎪ 01−− exp− ⋅ ⎬ ∏ 2⎨ 22nnσ πσ ⎪⎪ 0= ⎩⎭ 21N− 1xn⎧⎫ [] ⎪⎪ exp− ∏ 2⎨⎬ 22nn Mσ πσ = ⎪⎪ ⎩⎭ 0+ în care parametrul necunoscut este durata de propagare a semnalului de sondaj nτΔ= 00 44 22 Expresia plauzibilității poate fi adusă la o formă mai simplă 1N− 211⎧⎫ x ⋅ [] () ⎨⎬ ∑ 02/2;expNpn xn=− 22nσ 0= ⎩⎭ 2πσ () 01nM+− ⎧⎫ 21⎪⎪ +− [][][] () ⎨⎬ ∑ exp 2xnsn n s n n−−00 2− = ⎪⎪ 02nnσ ⎩⎭ în care primul factor nu depinde de parametrul necunoscut Pentru a găsi un estimator ML este necesar să maximizăm al doilea factor, ceeace echivalează cu minimizarea expresiei 01nM+− 2 −+ − = [][ ] [ ] () ∑ 002xnsn n s n n− 0nn= +− 0011nM nM+− 2 +− [][ ] [ ] ∑∑ 002xnsn n s n n−− 00nn nn== Al doilea termen este o constantă, egală cu energia semnalului de sondare − − 011nMM+ 22 −== [] [] ∑∑ 0ssnn smε 00nn m== expresiei analizate este echivalentă45Minimizarea cu maximizarea funcției − 01nM+ () [][ ] ∑ 00Jnxnsnn=− 0nn= Valoarea care maximizează funcția J este estimatorul ML a timpului ˆn 0ML Cu această estimare se determină estimatorul ML pentru distanța RADAR-țintă cΔ ˆ ˆRn= 0MLML 2 Estimarea parametrilor semnalului sinusoidal Vom relua problema estimării celor trei parametrii ai semnalului sinusoidal, al cărui model este 2 σ=+Φ+=−∼N [][][] () () 0cos 2 ; 0,1, , 1; 0,xn A fn wn n N wnπ Parametrul vector luat ]n considerare are forma T θ ⎡Φ⎤ 0Af= ⎣⎦ Cu datele x fixate, funcția de plauzibilitate, dependentă de vectorul θ este 1N− 211⎧ ⎫ ⎡⎤ x;θ exp cos 2pxnAfnπ =−−+Φ () [] () ⎨⎬ ∑ ⎦ 4602/22N⎣ = 02nσ ⎩⎭ 2πσ () 23 Pentru a maximiza plauzibilitatea este necesar să minimizăm, după θ, expresia 1N− 2 ⎡⎤ θ =Φ= − +Φ () [] ()() ∑ 00,, cos2JJAf xnA fnπ ⎣⎦ 0n= 1N− 2 =−Φ +Φ [] () ∑ cos cos 2 sin sin 2xnA fnA fnππ00 0n= Cu notațiile α=Φ=−Φ 12cos ; sinAAα forma funcției de minimizat, J, devine 1N− 2 απΦ= − − [] () () ∑ 01020,, cos2 sin2JAf xn fn fnαπ 0n= Sunt evidente relațiile − 222α Aarctgαα=+ Φ=; 12 α 471 Definim vectorii c (cosinus) și s (sinus) T ⎡⎤ π =⋅⋅−c () 00 01 cos 2 cos 2 2 cos 2 1ff fNππ ⎣⎦ T π=⋅⋅−s ⎡⎤ () 00 00 sin 2 sin 2 2 sin 2 1ff fNππ ⎣⎦ șiconstruim forma ⎡⎤ −⋅− ⋅ [] 12010xαα ⎢⎥ απ −− [] 10201cos2 sin2xffαπ ⎢⎥ ⎢⎥ α απ απ −− = − ⋅− ⋅xc s [] 12 1 0 2 02cos22sin22xf fα ⎢⎥ ⎢⎥ ⎢⎥ απ −− −− − () () [] 10 201cos2 1 sin2 1xN f N f Nαπ ⎣⎦ Prin calcul direct se verifică faptul că funcția J poate fi pusă sub forma Φ= ()() 0120,, , ,JAf J fαα T αα=− − − −xc sxc s ()() αα12 12 48 24 Facem notațiile pentru matricea H șivectorul α 10⎡⎤ ⎢⎥ 00cos 2 sin 2ffππ ⎢⎥ ⎡⎤ 1α ⎢⎥ ⋅⋅ === Hcsα [] ⎢⎥ ffππ00cos 2 2 sin 2 2 ; ⎢⎥ 2α ⎣⎦ ⎢⎥ ⎢⎥ −− () () 00cos 2 1 sin 2 1fN fNππ ⎣⎦ Cu acestea, funcția de minimizat, J, devine T =− −xHα xHα ()() () 120,,Jfαα şi apoi după frecvența digitală Am minimizat Vom minimiza mai întâi după α mai devreme o altă funcție 1T− J=θx-Hθ Cx-Hθ ()()() Dacă punem 1− CI u= stabilim, prin similitudine că 1− ˆTT= HH Hx () MLα Cu valorile determinate, rezultă pentru J o formă ce va fi minimizată după f T =− − xHαxHα ()() () 120ˆˆ ˆ ˆ,,Jfαα 11T−− TT TT⎛⎞ ⎛⎞ − =− x HHH Hx x HHH Hx () () ⎜⎟⎜⎟ ⎝⎠⎝⎠ 11−− ⎤⎡ ⎤ 49TTT TT⎡ =− xI HHH H I HHH Hx () () uu − ⎢⎥⎢ ⎥ ⎣⎦⎣ ⎦ Matricea factor din parantezele drepte este evident simetrică şi idempotentă 11−− TT TT⎡⎤ ⎡⎤ = IHHHH IHHHH () () uu−− ⎢⎥⎢⎥ ⎣⎦⎣⎦ 1111−−−− TT TT TTTT I HHH H HHH H HHH HHHH H () () () () u−−+ 1− TT IHHHH () u=− Cu acestea funcția de minimizat, J, ia forma 1− TTT⎡ ⎤ =− xI HHH Hx () () 120ˆˆ,,uJfαα ⎢⎥ ⎣⎦ 1− TT T T xx xHHH Hx =− () Minimizarea funcției J este echivalentă cu maximizarea termenului al doilea, primul fiind o constantă 1− TT T xHHH Hx () () 10Jf= 1− TT ⎛⎞ ⎤⎡⎤ cc T⎡ = xcs cs x ⎜⎟ [] [] ⎥⎢⎥ TT⎢ ⎜⎟ ss ⎢⎥⎢⎥ ⎣⎦⎣⎦ ⎝⎠ 1− TT T ⎤⎡ ⎤ cc cs c T⎡ = xcs x [] ⎢⎥⎢ ⎥ 50TT T sc ss s ⎢⎥⎢ ⎥ ⎣⎦⎣ ⎦ 25 Se stabilesc relațiile aproximative, dacă frecvența digitală nu este foarte aproape de 0 sau 0 5 11NN−− 1 TT π== = ≅cs sc ∑∑ 00 0cos 2 sin 2 sin 4 0fn fn fnππ 002nn== 11NN−− 21TN ==+≅cc () ∑∑ 00cos 2 1 cos 4fn fnππ 0022nn== 11NN−− 21TN ==−≅ss sin 2 1 cos 4fn fnππ () 00∑∑ 0022nn== 51 Prin calcul direct arătăm că funcția de maximizat este dublul periodogramei E suficient să determinăm frecvența la care apare maximul perodogramei pentru a avea estimatorul ML pentru frecvența digitală 1− N⎡⎤ 0 2T⎢⎥ xcs x [] () ⎢⎥ 10Jf= N 0⎢⎥ 2⎢⎥ ⎣⎦ T TT ⎛⎞ ⎡⎤ ⎡⎤ cc 2⎛⎞ xx = ⎜⎟⎜⎟ ⎢⎥ TT⎢⎥ ⎜⎟ N⎜⎟ ss ⎢⎥ ⎢⎥ ⎣⎦ ⎣⎦ ⎝⎠⎝⎠ T TT ⎡⎤ cx cx 2⎡⎤ = ⎢⎥ TT⎢⎥ Nsx sx ⎢⎥⎢⎥ ⎣⎦⎣⎦ 222 TT⎡⎤ =+ cx sx ()() N⎢⎥ ⎣⎦ 22 11NN−− ⎛⎞⎛⎞ 2⎡⎤ ⎢⎥ =+ [] [] ∑∑ ⎜⎟⎜⎟ 00cos 2 sin 2xn fn xn fnππ ⎢⎥ 00nnN== ⎝⎠⎝⎠ ⎣⎦ 2 1N− 22⎛⎞ − 022jfnπ xne= () ∑ ⎜⎟ 0Xf= 52[] nN0N= ⎝⎠ 26 21 If Xf= ()() 2Jf If= ()() N10 0 ˆ () 0arg maxMLfIf= 0f Frevența digitală la care apare maximul periodogramei calculată pentru datele 53 x din realizarea curentă,determină estimatorul ML pentru frecvență Având estimatorul pentru frecvență, se determină estimatorii ML pentru componentele vectoruluiα ⎤ 1T− cx 2TT⎡ ˆ= αHH Hx I = () uML⎢ T⎥ Nsx ⎢⎥ ⎣⎦ ˆα 1ML⎡⎤ ˆ= α ML⎢⎥ ˆα 2ML⎣⎦ 1N− ⎤ 2⎡ ˆ [] ∑ ⎢⎥ 0cos 2xnfnπ 0nN= ⎢⎥ = 1N− ⎥ 2⎢ ˆ [] ⎢⎥ ∑ 0sin 2xnfnπ 0nN= ⎣⎦ Se stabilește estimatorul ML pentru amplitudinea sinusoidei 1N− ˆ22 jfn− π 0222ˆˆ ===+ [] () ∑ 012ˆˆMLAxne Xfαα NN= 540n 27 În final, se stabilește estimatorul ML pentru faza inițială a sinusoidei 2ˆˆα − MLarctgΦ= 1ˆα 1N− ˆ [] ∑ 0cos 2xnfnπ 0n= =− 1Narctg− ˆ [] ∑ 0sin 2xnfnπ 0n= Determinarea direcţiei din care sosesc undele sonore în SONAR 55 Direcția din care vine un semnal sonor este legată de frecvențaspațială care are formula d == [] SfffHzβ00cos c Dacă sunt de estimat ML amplitudinea, frecvența digitală și faza inițială, suntem în situația studiată chiar mai înainte Vom presupune că frecvența analogică, geometria de dispunere a hidrofoanelor precum și viteza de propagare a sunetului în apă sunt cunoscute Pentru a estima frecvențaspațială, cu datele x construim periodograma 2 1N− 21jfnπ S− () [] ∑ SIf xne= 0nM= și determinăm, din poziția vârfului ei, estimatorul ˆ ,MLsf Se determină apoi estimatorul ML al unghiului β ˆ ,MLsfc ˆ arccosβ = 56ML 0fd 28